import geopandas as gpd
import trimesh
import numpy as np
import os
from shapely.ops import unary_union
from shapely.validation import make_valid, explain_validity
from shapely import is_valid
from shapely.errors import TopologicalError, GEOSException

def geojson_to_3d(geojson_path, output_path, scale=1000000, export_format='glb'):
    """
    将幸福经营所耕地格式的GeoJSON转换为3D模型
    网格命名规则：Id_姓名_LYJ
    适配CRS84坐标系统和包含Id、姓名、LYJ等属性的GeoJSON
    """
    # 检查输入文件是否存在
    if not os.path.exists(geojson_path):
        print(f"错误：输入文件不存在 - {geojson_path}")
        return

    # 创建输出目录
    output_dir = os.path.dirname(output_path)
    if not os.path.exists(output_dir):
        os.makedirs(output_dir)
        print(f"已创建输出目录：{output_dir}")

    # 验证导出格式
    supported_formats = ['glb', 'gltf', 'obj', 'stl', 'ply']
    if export_format not in supported_formats:
        print(f"错误：不支持的导出格式 {export_format}，使用默认格式 glb")
        export_format = 'glb'

    # 确保输出路径包含正确的扩展名
    if not output_path.lower().endswith(f'.{export_format}'):
        output_path = f"{output_path}.{export_format}"
        print(f"调整输出路径为：{output_path}")

    # 检查FCL是否可用（碰撞检测）
    try:
        import fcl
        fcl_available = True
        print("检测到python-fcl库，将使用精确碰撞检测")
    except ImportError:
        fcl_available = False
        print("未检测到python-fcl库，将使用边界框碰撞检测")

    # ---------------------- 1. 读取并解析GeoJSON ----------------------
    try:
        # 读取GeoJSON（适配CRS84坐标系统）
        gdf = gpd.read_file(geojson_path, crs="EPSG:4326")
        print(f"成功读取GeoJSON：共包含 {len(gdf)} 个地块，坐标系统为CRS84（EPSG:4326）")
    except Exception as e:
        print(f"读取GeoJSON失败：{str(e)}")
        return

    # 检查必要字段（核心命名字段）
    required_fields = ["Id", "姓名", "LYJ", "LCS", "DK_BM"]
    missing_fields = [f for f in required_fields if f not in gdf.columns]
    if missing_fields:
        print(f"错误：GeoJSON缺少必要字段 - {missing_fields}，无法按规则命名网格")
        return

    # 修复所有几何图形的拓扑错误
    print("开始修复几何图形的拓扑错误...")
    valid_geoms = []
    for idx, geom in enumerate(gdf.geometry):
        if geom is None or geom.is_empty:
            valid_geoms.append(None)
            continue

        try:
            if not is_valid(geom):
                # 打印无效原因（仅前10个）
                if idx < 10:
                    print(f"修复地块 {idx} 的拓扑错误: {explain_validity(geom)}")
                # 尝试修复几何
                fixed_geom = make_valid(geom)
                valid_geoms.append(fixed_geom)
            else:
                valid_geoms.append(geom)
        except Exception as e:
            print(f"地块 {idx} 修复失败: {str(e)}")
            valid_geoms.append(None)

    # 更新 geometries 并过滤无效几何
    gdf['geometry'] = valid_geoms
    gdf = gdf[gdf.geometry.notna() & ~gdf.geometry.is_empty]
    print(f"修复完成，保留 {len(gdf)} 个有效的几何图形")

    # 检查原始数据拓扑
    try:
        union = unary_union(gdf.geometry)
        if union.geom_type in ["MultiPolygon", "GeometryCollection"]:
            print(f"检测到 {len(union.geoms)} 个独立区域，可能存在原始数据重叠")
    except (TopologicalError, GEOSException) as e:
        print(f"警告：检查拓扑时遇到错误 - {str(e)}，跳过此检查")

    # 收集有效几何与坐标
    all_coords = []
    valid_geometries = []  # 存储有效地块的索引

    for idx, row in gdf.iterrows():
        # 1. 跳过空几何
        if row.geometry.is_empty or row.geometry is None:
            print(f"警告：Id为 {row['Id']} 的地块（姓名：{row['姓名']}）为空白几何，已跳过")
            continue

        # 2. 统一处理MultiPolygon/Polygon类型
        if row.geometry.geom_type == "MultiPolygon":
            polygons = list(row.geometry.geoms)
        elif row.geometry.geom_type == "Polygon":
            polygons = [row.geometry]
        else:
            print(f"警告：Id为 {row['Id']} 的地块（姓名：{row['姓名']}）几何类型为 {row.geometry.geom_type}，仅支持Polygon/MultiPolygon，已跳过")
            continue

        if not polygons:
            print(f"警告：Id为 {row['Id']} 的地块（姓名：{row['姓名']}）无多边形数据，已跳过")
            continue

        # 3. 检查核心属性是否为空
        core_attrs = {"Id": row["Id"], "姓名": row["姓名"], "LYJ": row["LYJ"]}
        empty_attrs = [k for k, v in core_attrs.items() if v is None or str(v).strip() == ""]
        if empty_attrs:
            print(f"警告：Id为 {row['Id']} 的地块缺少核心命名字段 {empty_attrs}，已跳过")
            continue

        # 4. 收集有效坐标（处理可能的Z坐标）
        has_valid_coords = False
        for poly in polygons:
            if poly.exterior is None:
                continue
            coords = list(poly.exterior.coords)
            # 移除Z坐标（如果存在），只保留X,Y
            coords_2d = [coord[:2] for coord in coords]
            if len(coords_2d) >= 4:  # 闭合多边形至少4个点
                all_coords.extend(coords_2d)
                has_valid_coords = True

        if has_valid_coords:
            valid_geometries.append(idx)
        else:
            # 构建临时名称用于日志
            temp_name = f"{row['Id']}_{row['姓名']}_{row['LYJ']}"
            print(f"警告：地块 {temp_name} 无有效顶点（需≥4个点），已跳过")

    # 检查是否有有效坐标
    if not all_coords:
        print("错误：未找到任何有效的坐标数据，无法创建3D模型")
        return

    # 计算坐标偏移量（CRS84为经纬度坐标）
    all_coords = np.array(all_coords, dtype=np.float64)
    min_lon, min_lat = all_coords.min(axis=0)[:2]  # 取经纬度最小值作为偏移原点
    offset = np.array([min_lon, min_lat], dtype=np.float64)
    print(f"坐标偏移原点（经纬度）：经度={min_lon:.8f}, 纬度={min_lat:.8f}")

    # ---------------------- 2. 构建3D网格（按“Id_姓名_LYJ”命名） ----------------------
    scene = trimesh.Scene()
    geometries = []  # 存储所有有效网格
    existing_meshes = []  # 用于碰撞检测的已创建网格

    # 配置三角化引擎
    try:
        import mapbox_earcut
        print("将使用 mapbox_earcut 作为三角化引擎")
        os.environ["TRIMESH_TRIANGULATE_ENGINE"] = "earcut"
    except ImportError:
        try:
            import triangle
            print("将使用 triangle 作为三角化引擎")
            os.environ["TRIMESH_TRIANGULATE_ENGINE"] = "triangle"
        except ImportError:
            print("警告：未找到第三方三角化引擎，使用trimesh默认引擎")
            os.environ.pop("TRIMESH_TRIANGULATE_ENGINE", None)

    # 核心参数
    EXTRUSION_HEIGHT = 1e-5  # 极小高度避免平面错误
    SIMPLIFY_TOLERANCE = scale * 0.00001  # 多边形简化阈值
    MIN_VERTICES = 3  # 简化后最小顶点数
    MIN_DISTANCE_THRESHOLD = 1e-4  # 碰撞检测最小距离

    for idx in valid_geometries:
        row = gdf.iloc[idx]
        # 提取核心属性（统一格式）
        Id = str(row["Id"]).strip()
        姓名 = str(row["姓名"]).strip().replace(" ", "_")
        LYJ = str(row["LYJ"]).strip()
        # 构建唯一网格名称
        base_mesh_name = f"{姓名}_{LYJ}"

        # 重新获取当前地块的多边形
        if row.geometry.geom_type == "MultiPolygon":
            polygons = list(row.geometry.geoms)
        else:
            polygons = [row.geometry]

        # 处理每个子多边形
        for sub_idx, poly in enumerate(polygons):
            # 子多边形名称（MultiPolygon时添加子序号）
            mesh_name = base_mesh_name if len(polygons) == 1 else f"{base_mesh_name}_子区{sub_idx}"

            # 检查多边形有效性
            if poly.exterior is None:
                print(f"警告：{mesh_name} 的子多边形无效（无外边界），已跳过")
                continue
            coords = list(poly.exterior.coords)
            coords_2d = [coord[:2] for coord in coords]  # 移除Z坐标
            if len(coords_2d) < 4:
                print(f"警告：{mesh_name} 的子多边形顶点不足（{len(coords_2d)}个），已跳过")
                continue

            # 修复无效多边形
            original_poly = poly
            try:
                if not is_valid(poly):
                    print(f"警告：{mesh_name} 的子多边形无效，尝试修复...")
                    poly = make_valid(poly)
                    # 修复后若为MultiPolygon，取面积最大的子多边形
                    if poly.geom_type == "MultiPolygon":
                        poly = max(poly.geoms, key=lambda p: p.area if p.exterior else 0)
                    # 重新获取修复后的坐标
                    coords = list(poly.exterior.coords) if poly.exterior else []
                    coords_2d = [coord[:2] for coord in coords]
                    if len(coords_2d) < 4:
                        print(f"警告：{mesh_name} 的子多边形修复后顶点不足，已跳过")
                        continue
            except Exception as e:
                print(f"警告：{mesh_name} 的子多边形修复失败：{str(e)}，使用原始多边形")
                poly = original_poly
                coords = list(poly.exterior.coords)
                coords_2d = [coord[:2] for coord in coords]
                if len(coords_2d) < 4:
                    print(f"警告：{mesh_name} 的子多边形顶点不足，已跳过")
                    continue

            # 坐标转换（经纬度→平面坐标）
            planar_coords = (np.array(coords_2d, dtype=np.float64)[:, :2] - offset) * scale

            # 简化多边形（仅复杂多边形）
            path = trimesh.load_path(planar_coords)
            original_vertex_count = len(planar_coords)
            if original_vertex_count > 500:
                simplified_path = path.simplify(tolerance=SIMPLIFY_TOLERANCE)
                simplified_count = len(simplified_path.vertices)
                if simplified_count >= MIN_VERTICES:
                    path = simplified_path
                    print(f"{mesh_name} 简化完成：顶点数 {original_vertex_count} → {simplified_count}")
                else:
                    print(f"{mesh_name} 简化后顶点不足，保留原始顶点")

            # 检查简化后的路径有效性
            if len(path.vertices) < MIN_VERTICES:
                print(f"警告：{mesh_name} 的子多边形顶点不足，已跳过")
                continue

            # 生成3D网格
            try:
                mesh = path.extrude(height=EXTRUSION_HEIGHT)
                if isinstance(mesh, list):
                    mesh = mesh[0] if mesh else None
            except Exception as e:
                print(f"警告：{mesh_name} 挤压生成网格失败：{str(e)}，尝试手动构建")
                # 手动构建平面网格
                vertices = np.column_stack([planar_coords, np.zeros(len(planar_coords))])
                faces = [[0, i, i+1] for i in range(1, len(planar_coords)-1)]
                if not faces:
                    print(f"警告：{mesh_name} 无法手动构建网格，已跳过")
                    continue
                mesh = trimesh.Trimesh(vertices=vertices, faces=np.array(faces))

            if mesh is None:
                print(f"警告：{mesh_name} 未生成有效网格，已跳过")
                continue

            # 碰撞检测
            collision_found = False
            if existing_meshes:
                current_bbox = mesh.bounds
                for existing_mesh in existing_meshes:
                    existing_bbox = existing_mesh.bounds
                    # 边界框快速检查
                    if not (current_bbox[0][0] < existing_bbox[1][0] and
                            current_bbox[1][0] > existing_bbox[0][0] and
                            current_bbox[0][1] < existing_bbox[1][1] and
                            current_bbox[1][1] > existing_bbox[0][1]):
                        continue

                    # 质心距离检查
                    current_centroid = mesh.centroid
                    existing_centroid = existing_mesh.centroid
                    distance = np.linalg.norm(current_centroid - existing_centroid)

                    if distance < MIN_DISTANCE_THRESHOLD:
                        if not fcl_available or (fcl_available and check_collision_fcl(mesh, existing_mesh)):
                            print(f"警告：{mesh_name} 与已有网格重叠，已微调位置")
                            mesh.apply_translation([1e-5, 1e-5, 0])
                            collision_found = True
                            break

            # 设置网格属性与名称
            mesh.metadata = {
                "name": mesh_name,
                "Id": Id,
                "姓名": 姓名,
                "LYJ": LYJ,
                "LCS": row.get("LCS", ""),
                "DK_BM": row.get("DK_BM", "")
            }
            mesh.name = mesh_name

            # 添加到场景与缓存
            geometries.append(mesh)
            existing_meshes.append(mesh)
            print(f"已创建3D网格：{mesh_name}（顶点数：{len(coords_2d)}）")

    # ---------------------- 3. 导出3D模型 ----------------------
    if not geometries:
        print("错误：未生成任何有效3D网格，无法导出模型")
        return

    try:
        scene.add_geometry(geometries)
        scene.export(file_type=export_format, file_obj=output_path)
        print(f"\n3D模型导出成功！路径：{output_path}")
        print(f"统计信息：共创建 {len(geometries)} 个网格，对应 {len(valid_geometries)} 个有效地块")
    except Exception as e:
        print(f"导出3D模型失败：{str(e)}")
        # 备选格式兜底
        fallback_format = "obj" if export_format == "glb" else "glb"
        fallback_path = output_path.rsplit(".", 1)[0] + f".{fallback_format}"
        try:
            scene.export(file_type=fallback_format, file_obj=fallback_path)
            print(f"备选格式导出成功：{fallback_path}（格式：{fallback_format}）")
        except Exception as fe:
            print(f"备选格式导出也失败：{str(fe)}")


def check_collision_fcl(mesh1, mesh2):
    """使用python-fcl进行精确碰撞检测"""
    try:
        import fcl
        # 构建FCL碰撞对象
        bvh1 = fcl.BVHModel()
        bvh1.beginModel(len(mesh1.vertices), len(mesh1.faces))
        bvh1.addSubModel(mesh1.vertices, mesh1.faces)
        bvh1.endModel()
        obj1 = fcl.CollisionObject(bvh1)

        bvh2 = fcl.BVHModel()
        bvh2.beginModel(len(mesh2.vertices), len(mesh2.faces))
        bvh2.addSubModel(mesh2.vertices, mesh2.faces)
        bvh2.endModel()
        obj2 = fcl.CollisionObject(bvh2)

        # 执行碰撞检测
        req = fcl.CollisionRequest()
        res = fcl.CollisionResult()
        fcl.collide(obj1, obj2, req, res)
        return res.is_collision
    except Exception as e:
        print(f"FCL碰撞检测出错：{str(e)}，默认判定为无碰撞")
        return False


# ---------------------- 执行转换 ----------------------
if __name__ == "__main__":
    # 输入：幸福经营所耕地格式的GeoJSON路径
    input_geojson = r"C:\Users\ma\Desktop\GIS2\幸福经营所耕地.geojson"
    # 输出：3D模型路径
    output_model = r"C:\Users\ma\Desktop\GIS2\幸福经营所耕地.glb"

    # 调用转换函数
    geojson_to_3d(
        geojson_path=input_geojson,
        output_path=output_model,
        scale=1000000,  # 缩放系数，保留经纬度精度
        export_format='glb'
    )
